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A GAS-KINETIC METHOD FOR HYPERBOLIC-ELLIPTIC EQUATIONS 
AND ITS APPLICATION IN TWO-PHASE FLUID FLOW 

KUN XU* 


Abstract. A gas-kinetic method for the hyperbolic-elliptic equations is presented in this paper. In the 
mixed type system, the co-existence and the phase transition between liquid and gas are described by the 
van der Waals-type equation of state (EOS). Due to the unstable mechanism for a fluid in the elliptic region, 
interface between the liquid and gas can be kept sharp through the condensation and evaporation process 
to remove the “averaged” numerical fluid away from the elliptic region, and the interface thickness depends 
on the numerical diffusion and stiffness of the phase change. A few examples are presented in this paper for 
both phase transition and multifluid interface problems. 
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1. Introduction. The study of the liquid-gas phase transition and the interface movement have both 
theoretical and practical interesting. The macroscopic governing equations for this phenomena are the mixed 
hyperbolic elliptic system, where the van der Waals-typc equation of state is usually used. Many numerical 
schemes have been proposed to solve the mixed type system, and the intensive investigation for the possible 
Riemann solver of the mixed type equations is still undergoing [24, 23, 6, 22, 5, 12, 14, 8]. Many criteria, such 
as viscosity capillarity and entropy rate admissibility conditions, have been well recognized in the capturing 
of realizable solutions. 

Physically, the van der Waals model can be rigorously derived from statistical mechanics, and the 
coexistence region of liquid and gas can be well predicted from the Maxwell construction. The particle 
interaction with nearby repulsion and long ranged attraction can naturally give the phase transition and 
surface tension properties [16, 7]. Based on the particle interaction pictures, many Lattice Boltzmann schemes 
have been developed, see [20, 21, 26, 9, 17, 3] and refences therein, and the particle interaction mechanism is 
used to capture both multifluid interface and phase transition process. Recently, combining the macroscopic 
van der Waals equation of state (EOS) and the mesoscopic Lattice Boltzmann method, He et. al. developed 
an interesting scheme for the capturing of liquid and gas interface and they successfully applied the scheme 
to the study of the Rayleigh- Taylor instability [10]. However, the scheme in [10] only applies the van der 
Waals EOS to the index function and it does not capture the liquid and gas phase change. Also, the densities 
of the liquid and gas in [10] are artificially assigned which may not be consistent with the values from the 
van der Waals EOS and the Maxwell construction. Based on different interface sharpening mechanism, such 
as the reinitialization in the level set method, many interface capturing schemes have been developed in the 
past years, see [25, 13, 15] and references therein. 

In this paper, we are going to develop a gas-kinetic BGK-type scheme for the hyperbolic elliptic system, 
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where the continuum and momentum equations are solved directly. The phase transition and the motion of 
multifluid interface are accurately captured by the current method. 


2. Governing Equations and Interface Capturing Mechanism. In the one- dimensional case, the 
governing equations for the isothermal hyperbolic-elliptic system are 


( 2 . 1 ) 


( <• ) + ( pu 

\pU h \pU 2 +p 


= 0 


where p and U are the density and velocity. For the multiphase flow and phase transition problems, the 
relation between the pressure p and the density proposed by van der Waals is quite satisfactory. The equation 
of state is 


pRT 2 

T=Tp- ap ' 


where R is the gas constant, T is the temperature, and a and b are constants. The critical temperature for 
the separation of liquid and gas is 


T c = 


8 a 

27 bR' 


When the fluid temperature is below the above critical value, phase segregation occurs. In this paper, we 
are going to study the fluids with fixed values a = 0.9, b = 0.25 and RT = 1.0. The corresponding critical 
temperature in this case is T c = 1.0666 /R. Since RT is less than T C R, the phase transition can appear in 
the current fluid system. The illustrative plot of the van der Waals EOS is shown in Fig. (5.1). The densities 
of liquid pi and gas p g can be obtained from Maxwell construction (equal area construction). The values in 
the plot arc l/pi = 0.494273, l/p g = 1.405065, l/p a = 0.574912, l/pp = 1.036251. The fluid density p can 
be catalogued in the following regions: 
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When the fluid density takes on values in the elliptical region, due to the negative slope of dp/ dp, 
fluid instabilities will be amplified. The fluid mixture in the elliptic region will evaporate to the gas state 
or condense to the liquid state. So, similar to the shock steepening mechanism, the van der Waals EOS 
has an intrinsic physical mechanism to separate different phases at the multifluid interface and sharpen the 
interface. This is the main reason that we can observe the sharp liquid gas interface in the real world. This 
property can also be used to develop an interface capturing scheme. 

Numerically, due to the cell size limitation and averaging process [28], the liquid and gas will be artificially 
mixed to form a numerical mixture with averaged density. If there is no any steepening mechanism at the 
interface, such as at the contact discontinuity wave of the compressible Euler equations, the thickness of 
the interface will grow with the square root of evolution time or total number of time steps. The real 
sharp interface between different fluids is enforced through the van der Waals-typc EOS or the molecular 
interactions. Once this kind of physics is incorporated in a multifluid numerical scheme, the numerically 
averaged density in the elliptical region will be condensed to the liquid or evaporated to the gas and the 
numerical interface can be kept sharp automatically. 
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In order to use the steepening mechanism at a multifluid interface, a scheme must be very accurate in 
predicting the liquid and gas densities first. In other words, even though there is no explicit terms about 
Maxwell equal area construction in Eq.(2.1), a scheme must have certain intrinsic dissipative mechanism, 
such as the implicit viscosity capillarity terms, and be able to pick up the physically correct solutions, such 
as the same liquid and gas densities from the Maxwell construction. Then, a numerically averaged density 
can be most likely located in the elliptic region. On the other hand, the numerical diffusion cannot be too 
large or it will overtake the physical steepening mechanism in the elliptic region. 

Currently, it seems difficult for any high-order scheme to predict a very accurate density jump at a 
multifluid interface. It is not surprising that many existing high-order schemes will give liquid and gas 
densities which probably depend on the interpolation limiters, CFL number, the cell size, even the Runge- 
Kutta time stepping techniques. In the current paper, we present a gas-kinetic scheme to solve Eq.(2.1). 
Due to the intrinsic diffusion and dissipative mechanism in the kinetic approach [28], the Maxwell equal area 
rule is implicitly achieved. The numerically obtained equilibrium densities of liquid and gas are very close to 
the theoretical values. At the same time, the phase boundary can be kept within two or three grid points. 
W ith this property, the kinetic method is used to simulate the evolution of multifluid interface, such as the 
merging of two liquid droplets. 

3. Gas-Kinetic Scheme for the Hyperbolic-Elliptic Equations. The BGK model has the stan- 
dard form 

(3i) u + uf, = 

r 

where / is the gas distribution function, g is the equilibrium state, and r is the particle collision time. Both 
/ and g are functions of space x, time and particle velocity u. 

In order to recover Eq.(2.1) from Eq.(3.1), the equilibrium state g can be constructed as 

g = Pi-)-e - A (“- c/ ) 2 , 

7 T 

where A is defined by 

A 

(3.2) 


and the variation of A is related to the density changes 
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where the functions A and D are well defined in the above equations. In the current paper, a fixed value 
RT = 1.0 is used. 

Due to the conservation properties in particle collision, / and g satisfy the compatibility condition 
(3-4) f (f - g)ip a du = 0 for ip Q = (1, u) T 

at any point in space and time. 
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The solution for the BGK model (3.1) is 


(3.5) 


f(x i+ x/ 2 ,t,u) = - f g(x',t',u)e 1 ),T dt' + e t/r fo{x j+ i/ 2 - ut), 

T Jo 


where x J+ \/ 2 is the cell interface and x' — Xj+\/ 2 —u(t — t') is the particle trajectory. There are two unknowns 
in the above equation. One is the initial gas distribution function f 0 at time t = 0, and the other is g in 
both space and time locally around (x j+1/2 ,t = 0). Similar to the BGK-type schemes for the Euler and 
Navier-Stokes equations [28], the numerical scheme based on Eq.(3.5), along with the compatibility condition 
(3.4), is described as follows: 

Step(l): Use MUSCL technique [27] to interpolate the conservative variables W = (p,pU) T , and obtain 
the reconstructed initial data inside each cell 


(3.6) 


Wj{x) = Wj{xj) + 


Wjfoj+i/a) ~ w j( x j- 1/2) 
x j+ 1/2 “ x j- 1/2 


(x — Xj) 


for x G [xj-\/ 2 ,Xj + i/ 2 ], 


where Wj(xj) is the cell averaged value, and Wj(xj_ 1 / 2 ) and Wj(xj + 1 / 2 ) are the values at the cell boundaries. 
A nonlinear limiter, such as van Leer’s limiter, is used in the current paper to get the cell boundary values. 
Prom the reconstructed data, the values p, V and their corresponding slopes, e.g. dp/dx and dU/dx, are 
known everywhere. Therefore, the variation A can be found subsequently through Eq.(3.3), such as dXjdx = 
D(p)dp/dx. 

Step(2): Based on the reconstructed data in Step(l), around each cell interface x j+i/2 , construct the initial 
gas distribution function /o, 


(3.7) 
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g l (l + a l (x — Xj + 1/2)) i x — x j+l/2 
g r (1 + a T (x - x j+1/2 )), x>x j+1/2 , 


where the states g 1 and g r arc the Maxwellian distribution functions defined in terms of the conservative 
variables at a cell interface. 


(3.8) 
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all coefficients in g l can be obtained as 
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Similar formulation can be found for g r . The coefficients a l r in Eq.(3.7) have the forms 


a lr = rn l { r + m 1 ^ u + m 2 r u 2 , 


which are derived from the Taylor expansion of a Maxwellian distribution function. The coefficients (m' 1 ’ r , mif , m' 3 ,r ) 
depend on (p‘, U 1 ), (p r ,U r ) and their corresponding slopes, i.e. 
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The detail relations are 
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Therefore, with the initially reconstructed data in Step (1), f 0 (x) in Eq.(3.7) are totally determined. For 
the sake of simplicity, we assume Xj+i/% — 0 in the rest of this paper. 

Step(3): The equilibrium state g is constructed as 


(3-12) g = g 0 (l + (1 - + H[x]a r j; -f At) , 

where H[ar] is the Heaviside function and g 0 is the state located at (x = 0, t = 0), 

(3.13) 9o = Po(-)*e- x o(u-u 0 )\ 

7 r 

The coefficients a / ,a r , and A in Eq.(3.12) have the forms 

a l,r = fn{ r 4- fn£ u -f m^ r u 2 , 

(3.14) A = A\ + A2U + A3U 2 , 


which have the same functional dependence on (dp/dx,dU/dx) and ( dp/dt,dU/dt ), as shown in Eq.(3.12). 

Taking both limits ( x — > 0) and (t — > 0) in Eq.(3.5) and (3.12), and applying the compatibility condition 
at (x = 0,t = 0), we can get macroscopic quantities Wo, 


( 3 - 15 ) w q ~ = J 9o^adu = J (^H[u] + g r ( 1 - H [u])) t/> Q du, 

where g l and g r are known from Step (2). Here Wq is the “averaged” flow variables at the cell interface, 
from which g 0 can be determined. Then, connecting W 0 to the cell centered values W J (x J ) and W j+1 (xj+i), 
we can get the slopes for mass and momentum distributions on both sides 
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for x < 0, 
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from which dp/dx.dU /dx.dX/dt can be obtained. Therefore, (a*,a r ) can be determined in a similar way 
as that in Eq.(3.12). The only unknown in Eq.(3.12) is A , which is related to dp 0 /dt , dU 0 /dt and d\ 0 /dt 
(= D(p 0 )dpo/dt) through the relations 
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To this point, we need to evaluate dpo/dt and d(pU)o/dt. 

Step(4): Substituting Eqs.(3.12) and (3.7) into the integral solution (3.5), we obtain the distribution func- 
tion / at x = 0, 

/( 0, f, u ) = jo 9o + Ji (a ! H[u] + a r (l - H[u])) ug 0 
(3.18) +J 2 M 0 + 73 ((1 - uta l )U[u]g l + (1 - uta r ){ 1 - H(u])p r ) , 


where 

70 = 1 - e- t/T , 

71 = r(— 1 + e~ t/T ) + te~ t/T , 

72 = T(tjr - 1 + e” t/,r ), 



The only unknown in Eq.(3.18) is A , which is a function of (dpo/dt, dU 0 /dt), see Eq.(3.18). Since the 
compatibility condition must be satisfied everywhere in space and time, it can be integrated in a whole CFL 
time step AT at x = 0 


(3.19) 



(/(0, t, u) - g( 0, t, u)) ipadtdu = 0, 


from which 


(3.20) 
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= J [-r 3 So + Tiu (a'H[u] + a r (l - H[u])) g 0 
+ r 3 (HMs' + ti-HH)^) 

+ T 4 u (a'H[w] 3 ( + a r (l - H[u])p r )] ip a du 


arc obtained. All terms on the right hand side of the above equation are known, and 

T 0 = AT - r(l - e~ AT/T ), 

Ti = t (-AT + 2t(1 - e' AT/r ) - A Te at/t ) , 

r 2 = ^AT 2 - tAT + r 2 ( 1 - e~ AT/r ), 
r 3 = r(l-e- AT/r ), 

r 4 = -r (-A Te' AT/r + t(1 - e^ AT/r )) , 

r 5 = r(AT-r(l-e- AT/T ))- 

Thus, (dpo/dt. dUo/dt), as wel as d\ 0 /dt, can be obtained from Eq.(3.21). 

Step(5): The time-dependent numerical fluxes of mass and momentum across the cell interface is 


(3.21) 

where / is given in 


Fwj+ 1/2 — ( j? P — f u V ; a/j+i/ 2 ( 0 , t, u)du, 

Eq.(3.18). The update of flow variables inside each cell becomes 

1 r At 

W7 +1 = Wj 1 + — J (Fwj-I / 2 - F\Vj+i/2)dt- 
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4. Numerical Examples. In this section, we arc going to present a few test cases in both ID and 2D. 
The van Leer limiter is used for interpolations of p and pU at the beginning of each time step. The time step 
AT is determined by the Courant- Friedrichs- Levy condition with CFL number =0.25. The collision time r 
is defined as 

(4.1) r = CiAr + C 2 Ar ^l~ P ^ , 

pi 

where p l = A (p l ), p r = A (p r ), and C x = 0.05, C 2 = 2.0 are fixed in all calculations. 

4.1. ID Shock Tube Cases. In the following, four shock tube cases presented in Shu’s paper [22] are 
tested using the current kinetic scheme. 

CASE(l) 

The initial condition for this case is the exact liquid and gas densities from the Maxwell construction, 

{1/PL = 0.494273, U L = 1.0)| x<0 . 5 and (l/p R = 1.405065, U R = 1.0)| x > 0 . 5 . 

The cell size used is Ax = 1/200. This test case is mainly to see if the scheme can keep the admissible jump 
from the Maxwell construction. The numerical result (solid line and circles) for the density distribution is 
shown in Fig. (5.2), where the dashed and dotted lines represent 1/p/, 1 j p g , 1 / p Q and l/pp respectively. 

CASE(2) 

The second case has the following initial condition, 

(1 1 Pl = 0.54, Ui = 1.0) and (1 fp R ^ 1.8517, U R = 1.0). 

This initial jump satisfies the Rankine-Hugoniot condition, but does not satisfy the physical principles, 
such as viscosity capillarity condition. With the cell size Ax = 1/200, our simulation results are shown in 
Fig. (5.3), where the solid line is obtained with a much refined mesh Ax = 1/2000. This case clearly shows 
that the current scheme can pick up the physically admissible solution. There is no oscillations at the liquid 
phase around the interface. Our results arc favorable in comparison with the ENO-typc method [22]. 
CASE(3) 

This case has the initial condition 

(1/pL = 0.45, Ul = 1.0) and (1 fp R = 2.0, Ur — 2.0). 

Fig. (5.4) shows the results with cell sizes Ax = 1/200 and Ax = 1/2000. From this figure, we can also 
observe the sharp interface between the liquid and gas phases. 

CASE(4) 

The initial condition for this case is 

(1/p, U) = (0.8 + 0.2sin(x), 1 — 0.5cos(x)). 

The initial density is entirely in the elliptic region. Periodic boundary condition is used. The solutions 
with cell size Ax = 1/400 at different output time are shown in Fig. (5. 5). These figures clearly show the 
flow instability in the elliptic region and how the densities eventually go to the well defined liquid and gas 
densities, even though the Maxwell construction is not explicitly used in the current scheme. For the liquid 
and gas phases, the numerical densities obtained are about 1/p = (0.49400, 1.40175). The differences between 
the numerics and the theoretical values (0.494273, 1.405065) are less than 0.5%. This is a very good case to 
test the ability of any high-order scheme to capture the correct density jumps around the phase boundary, 
as well as the sharpness of the interface. Our scheme can capture the jump within 2 or 3 cells, as shown in 
Fig(5.5d). 
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4.2. Liquid-Gas Interfaces in 2D Cases. In 2D cases, in order to capture the movement of a mul- 
tifluid interface the inclusion of surface tension and gravity becomes important. In the test case (5), the 
gravitational force pG is implemented in the y-momentum equation for the liquid phase, and the nondimen- 
sional magnitude of G is assigned the value 0.25. In the test case (6), an additional body force /cpVV 2 p 
is added in the momentum equations to recover the surface tension effect [19, 3, 10]. The nondimensional 
coefficient n used in case (6) is equal to 5.0 x 10 . 

CASE(5) 

This is a dam break problem. Many schemes have been used in this kind of free surface problems. 
A short list includes Volume of Fluid (VOF), Boundary Integral Techniques, Front Tracking Method, and 
Arbitrary Lagrangian-Eulcrian (ALE) method, see [11, 29, 4, 1] and references therein. The cell size used 
in our study is Ax = A y = 1/100. The schematic construction for this problem is shown in Fig. (5.6), 
where the densities of the liquid and gas are assigned with the values from the Maxwell construction, i.e. 
\/pi = 0.494273 and 1 j p g = 1.405065. The initial velocity of both gas and liquid arc zero, and no surface 
tension is included in this case. Due to the numerical diffusion, any index function to describe the liquid 
and gas interface will get smeared in the Eulerian advection scheme, and the smearing is proportional to the 
square root of the number of time steps used. This is basically the main reason for the level set method to use 
reinitialization to keep the interface sharp [2]. However, in our case, since the van der Waals EOS is used to 
describe the liquid and gas phases, any smeared density at the interface is most likely to locate in the elliptic 
region and the flow instability in these region will automatically steepen the interface. More specifically, 
the condensation and evaporation process around phase boundary could move the averaged density to the 
liquid or gas phases, and this effect compensates the numerical dissipation in the advection scheme. Fig. (5.7) 
shows the time evolution of the liquid-gas interface, and the interface thickness keeps two or three mesh size 
regardless the time steps used to get the final results. Fig. (5.8) shows the locations of the leading liquid 
front. The numerical results arc compared with the experimental data in [18]. From this figure, we observe 
that the numerical speed is slower than the experimental speed. The reason for the difference is that in the 
current calculation the density ratio between liquid and gas is about 2.8, and the experimental data was 
obtained for the water and air, and their density ratio is about 800. Therefore, the relative aerodynamical 
resistencc is much higher in the current study. Fig. (5. 7) shows a very interesting phenomena that there is 
a bore shock and a rarefaction wave in the liquid phase. This solution is amazingly close to the solution by 
solving the shallow water equations. In other words, the current direct numerical simulation in some sense 
validates the approximation used in theoretical derivation of shallow water equations. 

CASE(6) 

This test case is about the collision of two droplets. Similar to the last case, the initial densities of 
the liquid and gas phases are assigned with the theoretical values again from the Maxwell construction, i.e. 
1/p, = 0.494273 and 1 / p g — 1.405065. The cell size used in this case is Ax = Ay = 1/100. The initial 
droplets with radius R = 0.055 are moving toward each other with a velocity magnitude of U = 0.125. 
No gravity is included in this case. Fig. (5.9) shows the time evolution of the droplets. The collision and 
merging of the droplets can be observed. Due to the steepening mechanism at the fluid interfaces from 
the van der Waals EOS, the sharp interface is kept in the time evolution process. Fig. (5. 10) shows the 
density distribution across the central lines in both the x- and the y-directions of Fig. (5.9i) , where the phase 
boundaries keep 2 mesh points even though 1600 time steps have passed at that output time. 

5. Discussion and Conclusion. In this paper, we have developed a gas-kinetic scheme for the 
hyperbolic-elliptic system, where the van der Waals equation of state is used to describe the phase transition. 
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At the same time, the instability or sharpening mechanism in the elliptical region is used to capture the fluid 
interface. Many test cases validate the current approach. 

The current paper is only a first step in the study of multifluid flow by solving the mixed type system. 
Hopefully, in the near future we can answer the following questions. 

1. If we only need to describe the motion of multifluid interfaces without including phase transition, such 
as problems related to incompressible multifluid, we need to find some ways to simplify the van der Waals 
EOS in order to get a simple (not simpler) numerical method. But, the property dp/dp < 0 has to be kept 
so as to get the interface sharpening mechanism. 

2. The current approach can be used to capture the interface breaking and merging. Similar to the level 
set method, the interface topological changes is easily handled. However, since the equations are basically 
describing the phase transition problems, the mass conservation for each individual phase cannot be guar- 
anteed in the evolution process. How to improve the mass conservation property for each individual phase 
under the current framework is an important question that needs to be addressed. 

3. In order to increase the density ratio between the liquid and gas phases, we need to use a more realistic 
EOS to cope with that. 

The preliminary results presented in this paper arc very promising and encouraging. Since there is no 
tracking, index function, or special treatment used around the multifluid interfaces, the extension of the 
current method to simulate hundreds or even thousands of bubbles (or droplets) becomes possible. Also, 
the implementation of interface physics in the capturing of interface movement should be a reasonable and 
reliable approach. This kind of scheme may provide a new way to solve many multifluid engineering problems 
after its further developement. 
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FiG. 5.4. Distribution of l/p. Circles are the simulation result obtained with cell size Ax = 1/200. The solid line is the 
result obtained with a much refined mesh Ax — 1/2000. 



FiG. 5.5. The solid lines are the distributions of 1/p at different output times. The mesh size used is Ax = 1/400. (a) 
t=0. 0, (b) t=0.1 } (c) t=1.0, (d) t= 100.0. Circles are added in the plot (d) to show the number of grid points around the 
multifluid interfaces. 
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Fig. 5.7. Lxquid-gas interfaces at different output time, (a) t = 0.0, (b) G/a = 0.5, (c) tyj G/a = 1.0, (d) tyj G/a = 



Fig. 5.8. The horizontal axis is tyjG/a and the vertical axis is x/a, where x is the location of leading liquid front. The 
solid line is the time evolution of the leading liquid front. The density ratio between liquid and gas is around 2.8. The circle 
is the experimental data in [ 18 ], where real water and air with density ratio around 800 were used. 
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Fig. 5.9. Time e^oiufion of the collision of two droplets. The output times are (a) t=0, (b) t=0.2, (c) t—0.25, (d) t=0.3, 
(e) t= 0 . 40 , (f) t—0. 60, (g) t=0.80, (h) t-1.20, (i) 1=1.60. 
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FlG. 5.10. The distribution X/p. (a), along the central line of Fig. (5. 9i) in the x-direction, (b). along the central line 
of Fig.(5.9i) in the y-direction. Since both the liquid and gas are treated as the compressible flow in the current study, small 
density fluctuations appear in the dynamical transport process , especially in the gas phase. 
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